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Abstract 



The fundamental problem of optical wave propagation is the determination of the 
field at an observation point, given a disturbance specified over some finite aper- 
ture. In both vacuum and inhomogeneous media, the solution of this problem is 
given approximately by the superposition integral, which is a mathematical ex- 
pression of the extended Huygens-Presnel principle. 

In doing so, it is important to find the atmospheric impulse response (Green's 
function). Within a limited but useful region of validity, a satisfactory optical 
propagation theory for the earth's clear turbulent atmosphere could be developed 
by using Rytov's method to approximate the Helmholtz equation. In particular, 
we deal with two optical problems which are the time reversal and apodization 
problems. The background and consequences of these results for optical commu- 
nication through the atmosphere are briefly discussed. 
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Chapter 1 



Introduction 



In contemporary radiophysics, atmospheric optics and wireless communication, 
one often studies the propagation of electromagnetic waves in the atmosphere. In 
doing so, it is increasingly important to take into account the turbulent state of 
the atmosphere, a state which produces fluctuations in the refractive index of the 
air. 

Determining the turbulence effect on wave propagation requires one to solve Maxwell's 
equations with a refractive index that is stochastic. To solve for effects of turbu- 
lence in the low-power regime, one starts with Maxwell's equations 





(1.1) 



dt 




(1.2) 
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deir t) 

V X B{f, t) = HoPir, t) + Hoe{r, ^)— (^rid (1.3) 

V-5(r,t) = 0, (1.4) 



where e(f, t) is the electric field, S(r, t) is the magnetic field, e(r, t) is the permit- 
tivity, jjjQ is the magnetic permeability, p{f, t) is the current density, and g(r, t) is 
the charge density. Assume there is a time-harmonic variation of the electric field 

e{r,t) = exp{-icut)E{r) (1.5) 

where E{r) can also vary with time but the time scale of variation is much longer 
than the harmonic variation. The air density affects the refractive index as 

e(r) = eon2(r) (1.6) 

Then the wave equation describing the electric field derived from the above rela- 
tions is 



AE(r) + ky(r)E(r) - 2ih^^!fi^E^ + 2V{E(r) ■ Vln[n(m = (1.7) 

c ot 



where ko — — ^ is the free space wavenumber with c being the speed of fight 
in vacuum. 
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By reasonable simplifications (details in Chapter 4) , one finally need to study the 
Helmholtz equation (reduced wave equation) for the turbulent medium 



AE{r)+kln'^{r)E{r) =0 



(1.8) 



with appropriate boundary conditions describing the evolution of the harmonic 
amplitude for the components of the electric field. 

Within a limited but useful region of validity, a satisfactory optical propagation 
theory has been developed by using Rytov's method to approximate the Helmholtz 
equation. If it is assumed that the magnitude of the air-density inhomogeneity is 
small, the refractive index is then 



where 5ni(r) << 1. 

By using the Rytov approximation, we can obtain not only the solution, but also 
the Green's function of the Helmholtz equation for the turbulent medium with the 
aid of the Green's function for free-space propagation. 

It can also be shown that the Green's function in turbulence obtained by Rytov's 
method satisfies the same reciprocity condition as the free-space Green's function. 
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A better understanding of how randomness and turbulence affect the propagating 
wave fields can lead to better solutions to many problems in the optical society. 
Two examples that we are particularly interested when understanding wave tur- 
bulence interaction is important are: 

I. ) Time reversal problem. 

II. ) Apodization problem. 

In both these two examples we are concerned with the relationship between the ini- 
tial distribution of light over the exit pupil of an optical system and the amplitude 
distribution in the object plane. This is based on the extended Huygens-Fresnel 
Principle in random medium developed by H.T.Yura [12, 13]. The superposition 
integral requires an explicit form of the atmospheric impulse response (Green's 
function), which we have obtained by using Rytov's method. Further analysis is 
given for each problem. 

The structure of this thesis is as the following: 

We start Chapter 2 by constructing the Green's function for the free-space 3-D 
Helmholtz equation ()2.3p with the Sommerfeld radiation condition ()2.4|) . With a 
point source at a point ^, the Green's function is represented by ()2.15p . 
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In Chapter 3 we will give a brief introduction to the refractive index variation. We 
will also explain how the refractive index variation is determined by the tempera- 
ture fluctuations and present some experimental results. 

In Chapter 4 and 5, we first introduce the Rytov approximation method, and then 
derive the solution and Green's function for the Helmholtz equation in turbulent 
medium step by step. At the end of Chapter 5, the reciprocity property of the 
Green's function for general Helmholtz equation will be proved. 

Further simplification by parabolic approximation is conducted in Chapter 6. 
When solving for the solution of Helmholtz equation, we specialize the initial con- 
dition into two cases: plane wave and beam wave. An explicit form of the Green's 
function will also be derived by parabolic approximation. 

In Chapter 7 and 8, we will apply the results obtained in previous work to two 
important applications: time reversal problem and apodization problem. 



A brief discussion is given in Chapter 9. 
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Chapter 2 

Green's Function for Free-Space 
3-D Helmholtz Equation 

We consider the 3-D Wave Equation 

-{Uu-c^AU)=q{f,t) (2.1) 

where q{f,t) is a source. If q{f,t) = q{f)exp{—iujt) represents a source oscillating 
with a single frequency a;, then the entire motion reduces to a wave motion with 
same frequency uj after an initial transient period, so that we can write 



U{f,t) = U{r)exp{—iu>t) 



(2.2) 
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Thus ()2.1j) reduces to the 3-D Helmholtz equation 

-iA + k')Uir) = f{r} (2.3) 

where k = ^, f{r) = c^^q{r) 

The function U{r) satisfies this equation in some domain D <Z R with boundary 
dD^ and it also satisfies some prescribed boundary conditions. We also assume 
that U (r) satisfies the Sommerfeld radiation condition 

lim riUr - ikU) = (2.4) 

r— >oo 

which simply states that the solution behaves like outgoing waves generated by 
the source. 

We construct a Green's function G{r,^) for (|2.3|) . so that G{f,^) satisfies the 
equation 

-iA + k^)G = S{x)6{y)6{z) (2.5) 



Using the spherical polar coordinates, the 3-D Laplacian can be expressed in terms 
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of radial coordinate r, 

AG = Grr + -Gr (2.6) 

r 

so that ()2.5|) assumes the form 



1 rl riG 

-[-|-(^'?^) + ^'^]='^W, 0<r<oo (2.7) 



with the radiation condition ()2.4|) . 



For r>0, the function G satisfies the homogeneous equation 



4^(r^^) + ^^G = (2.^ 



Or, equivalently, 



o2 

— (rG) + A:2(rG) = (2.9) 



This equation admits a solution of the form 



rG(r) = Ae'"' + Ee-*"^" (2.10) 



or 



gjfcr ^—ikr 



G{r) = A— + B—— (2.11) 



where A and B are arbitrary constants. In order to satisfy the radiation condition, 
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we need to set B=0, and hence, solution ()2.11|) becomes 



pikr 

Q(r) = A (2.12) 

r 



To determine A, we use the spherical surface of radius e and divergence theorem, 
so that 

r BG r A 1 

lim / —dS = lim / -e'^''(ik - -)dS = -1 (2.13) 



from which we find A = ^ as e ^ 0. Consequently, the Green's function takes 
the form 

pikr 

= (2-14) 

Physically, this represents outgoing spherical waves radiating away from the source 
at the origin. With a point source at a point ^, the Green's function is represented 
by 

exp{.fc|f-e1} ^2.15) 
^ ^ 47r|f-e| 

where r and ^ are position vectors in R^. 



Finally, when k=0, this result reduces exactly to the Green's function for the 3-D 
Poisson equation. 
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Chapter 3 

Introduction to Refractive Index 
Variations 

We are all familiar with the twinkling of stars at night and the shimmering of dis- 
tant objects on a hot day. These visual effects are caused by small refractive index 
inhomogeneities in the atmosphere, which in turn are produced by atmospheric 
turbulence. Here we refer the word turbulence to the density fluctuations arising 
from atmospheric temperature fluctuations. When a beam of light passes through 
the air above a road heated by the sun, the randomly fluctuating air temperature 
produces small refractive index inhomogeneities that affect the beam propagation. 

Consider for example an initially well-deflned phase front propagating through a 
region of atmospheric turbulence. Because of random fluctuations in phase veloc- 
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ity the initially well defined phase front will become distorted. This alters and 
re-directs the flow of energy in the beam. As the distorted phase front progresses, 
random changes in beam direction and intensity fluctuations occur. The beam is 
also found to spread in size beyond the dimensions predicted by diffraction theory 
in free space. 

The reason why all this happens, as we have stated, is atmospheric turbulence 
that arises when air parcels of different temperatures are mixed by wind and con- 
vection. The individual air parcels, or turbulence cells, break up into smaller cells 
and eventually lose their identity. In the meantime, however, the mixing produces 
fluctuations in the density and therefore in the refractive index of air. A cartoon 
of laser beam propagation in turbulence is given in Figure ()3.1|) 

To describe these random processes, we need to find a way to define the fluctuations 
that are characteristic of turbulence. This is done by expressing each quantity as 
the sum of an average and a fluctuation term. We can write the random variation 
of the refractive index as 




with no =< n{r) > 



(3.1) 
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TURBULENCE CELLS (Large) 
(a) 



SCREEN 



BEAM 




TURBLILENCE CELLS (Small) 
(b) 



SCREEN 



Figure 3.1: (a) Laser beam is deflected by turbulence cells that are larger than the beam 
diameter, and (b) Laser beam is broken up by turbulence cells that are smaller than the 
beam diameter. 



Here the sharp brackets indicate the ensemble average 



< n(r) >= - r n(f,t)dt 

T Jo 



(3.2) 



where r is a time that is large compared to the lifetime of a fluctuation. A rep- 
etition of the averaging operation has no effect on < n{r) >, and it follows that 
< [< n{r) >] >=< n{r) > and < 5n{r) >= 0. 



In later chapters when we consider the Helmholtz equation in turbulence (details 
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are given in Chapter 4) 

Au{r) + kln^{r)u{r) = (3.3) 

We can write 

n{r) = 1 + Sni{r) (3.4) 

since we can assume < n{r} >= 1 by absorbing a constant term into the wavenum- 
ber ko, and write 

5n{r) = S ■ ni{r) (3.5) 
where 6 ■ ni{r) << 1, and 6 was inserted into this expression to show smallness. 

Now let's return to the original expression 1)3.11] . The refractive index fluctuation 
Sn{r) is proportional to the density fluctuation 6N{r), which in turn is governed 
by temperature fluctuations. The dependence of 5n{'f) on 6N{f) follows from the 
Lorenz-Lorentz law 

where Ci is a constant (for a given frequency of the electromagnetic wave) and 
is the particle density. For air, with n ~ 1, equation ()3.6|) may be approximated 
by 

(n - 1) = ^NC, (3.7) 
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Differentiation of the above equation yields tlie relationsliip between 5n and 5N: 



If we assume tliat air obeys tlie perfect gas law, then N = and it follows from 
(K7\i that 



where p is pressure expressed in millibars and T is temperature measured in 
Kelvins. For the optical region of the spectrum, the constant factor is given 
approximately by 79 x 10^®, thus 



We see that the refractive index depends on both pressure and temperature. Since 
pressure variations are relatively small, fluctuations of the refraction index are 
primarily due to temperature fluctuations. Because their short lifetime and small- 
scale sizes may be treated as an adiabatic process, hence 






(3.10) 



J. 



(3.11) 



pT(i-T) = constant 



where 7 ^ 1.4 for air. 



CHAPTER 3. Introduction to Refractive Index Variations 15 
Differentiation of equation ()3.11j) leads to 

5p , 7 ,(5T 



p 7 — 1 T 



(3.12) 



Now witli tfie aid of ()3.10|) we can eliminate 6p and find that the refraction index 
fluctuations are proportional to the temperature fluctuations, that is 



79 X 10^^ V 

6n = ^,6T (3.13) 

7 — 1 i ^ 



The single most important parameter appearing in almost all equations that de- 
scribe beam disturbances caused by turbulence is the refractive index structure 
coefficient C„. It is governed by the pressure and temperature difference at two 
points separated by the distance r (measured in centimeters) and is given by 



Cn = [79 X 10-'^]Ct (3.14) 



where the temperature structure parameter is 



Ct = V< (^1 - ^2)2 > r-3 (3.15) 

The temperatures T, Ti, and T2 are all in Kelvins, and p is the atmospheric pressure 
in millibars. For strong turbulence, C„ = 5 x 10~^; for intermediate turbulence, 
C„ = 4 X 10^®; and for weak turbulence, C„ = 8 x 10~^. These are typical values. 
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The structure parameter usually appears in equations in the form C^, which varies 
from about 10~^^m~3 for extremely weak turbulence to 10~^^m~^ or more when 
the turbulence is strong. This latter value is usually observed near the ground 
in direct sunlight. Measurements of C„ with temperature sensors 1.6m above 
the ground have shown that the minimum value of C„ occurs about one to two 
hours before sunrise and after sunset. The peak values have been recorded around 
noon on sunny days. As might be expected form equation ()3.14j) . the structure 
parameter C„ decreases with altitude. This is shown in Table 3.1. 

Table 3.1 Typical values of C„ as a function of height. 



Height (Km) 


C„(m-i/3) X 10^ 


0.001 


30 


0.003 


20 


0.01 


15 


0.03 


10 


0.1 


6 


0.3 


4 


1.0 


1 


3.0 


1 



The movement of small index-of-refraction inhomogeneities through the path of 
beam causes random deflection and interference between different portions of the 
wavefront, which can lead to an internal breaking up of the beam spot into smaller 
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hot spots. Figure ()3.2p shows a typical instantaneous intensity distribution of a 
He-Cd laser beam at 15.5 km. The initial intensity distribution of this beam at the 
laser was Gaussian. The photo illustrates clearly the complexity of the transmitted 
pattern that consists of amorphous areas of random sizes and shapes. Other pat- 
terns may have a reticulated appearance with sharp lines bounding large polygonal 
areas. The dominant size of the hot spots in the pattern is given approximately 
by V LA, where L is the distance from the laser to the observation screen. These 
bright patches of about 1-cm diameter are typical for L ^ 1km. For space-to-earth 
paths, L would be the distance between the observation screen and the turbulence 
layer in the atmosphere. 




Figure 3.2: Instantaneous intensity distribution of a laser beam at 15.5 km. Area of 
screen covered by photo is about 1.20 m by 1.55 m. 
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0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4 1.6 1.8 2.0 2.2 



TIME (Seconds) 

Figure 3.3: Typical received signal fluctuations for a small detector 145 km from the 
transmitter. (From A. L. Buck, Applied Optics 6, 703, April 1967) 

Because of the constantly changing pattern, a small detector placed in the beam 
will measure intensity fluctuations or scintillation. Typical power fluctuations 
recorded by a detector whose aperture size is considerably less than the beam 
diameter are shown in Figure ()3.3|1 . 

The temporal frequency of the intensity fluctuations recorded at a fixed point 
within the beam usually varies between 1 and 100 Hz. 

In determining the turbulence effect on wave propagation, it was Rytov who pro- 
posed an approximation that included diffraction effects. This resulted in a wide 
range of validity and many problems have been solved as interest in wave propaga- 
tion in turbulence expanded because of new technologies and areas of application. 
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These areas include sending laser beams efficiently through the atmosphere, re- 
moting sensing of the atmosphere on earth and correcting atmospheric distortion 
to allow better resolution in optical images. 

When propagating laser beams parallel to the ground, one found that scintillation 
increased with increasing turbulence as the Rytov theory predicted until a certain 
level was reached at which the measured scintillation saturated. Fortunately, even 
when the scintillation is saturated Rytov's method typically gives the correct an- 
swer for phase disturbances, thus allowing one to treat many problems of practical 
interest. It has generally been found that the Rytov approximation gives a good 
approximation if 



where ko is the free-space wavenumber, L is the distance from the laser to the 
observation screen and C„ is the refractive index structure coefficient. 




(3.16) 
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Chapter 4 

Rytov Approximation I - Solution 
of Helmholtz Equation 

Determining the turbulence effect on wave propagation requires one to solve Maxwell's 
equations with a refractive index that is stochastic. 

Assume there is a time-harmonic variation of the electric field 

e{f,t) = exp{-iujt)E{r) (4.1) 

— * 

where E{r) can also vary with time but the time scale of variation is much longer 
than the harmonic variation. 

The wave equation describing the electric field derived from the Maxwell's equation 
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is then 



AE{r) + kln^{r)E{f) - 2^— ^^""'^^^^""^^ + 2V{E{r) ■ Vln[n{f}]} = (4.2) 



where fco = ^ = ^ is the free space wavenumber with c being the speed of hght 



in vacuum, and n(r) is the refractive index. 



If the rate of change of E{r) with time is much less than the sinusoidal variation 
exp{—iut) , then the third term is negligible compared to the first two. For optical 
frequencies this condition is satisfied. 



In addition, if the propagation wavelength A is much less than the inner scale, 
we can ignore depolarization effects that is due to the last term in (j4.2j) . This 
corresponds to scatterers being large relative to the wave length and the direction 
of the scattered wave very close to the direction of the original scattered wave which 
is the case for the regime that we will consider. For visible wavelengths, typical 
inner-scale sizes of a millimeter satisfy this condition, in which case one obtains a 
scalar equation for each of the electric field components separately. Therefore in 
this regime the Helmholtz equation is 



AE{f) + kln\r)E{f) =0 



(4.3) 
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with appropriate boundary conditions describing the evolution of the harmonic 
amplitude for the components of the electric field. 



The equation for one component is 



Au{f) + k^n\f)u{f) = (4.4) 



In the Rytov method the solution is expressed as 



u{r) = expmf)] (4.5) 



This leads to the nonlinear Riccati equation 



A$(r) + |V*(r)|' = (4-6) 



If it is assumed that the magnitude of the air-density inhomogeneity is small, the 
refractive index is then 

n{f) = 1 + Sni{r) (4.7) 
where (5ni(r) << 1, and 5 was inserted into this expression to show smallness. 



By perturbation method, we assume that the solution can be written as a power 
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series in 6, i.e. 

$(f, 5) = ^o{f) + 5^i{f) + 52$2(f) + (4.8) 

By inserting ()4.7|) and ()4.8|) into ()4.6|) . and by separating the equations based on 
the power of 6, one obtains a system of equations 

6^ : A$o(f) + V$o(r) ■ V$o(r) = -kl (4.9) 

5^ : A$i(r) + 2V$o(f) ■ V$i(f) = -2A;^ni(f) (4.10) 

: A<l>2(r) + 2V$o(r) ■ V$2(r) = -A;^ni(r) - V$i(f) ■ V$i(r) (4.11) 



m— 1 

6"^ : A<l>™(r) + 2V<l'o(r^ • V<l'„(f) = - ^ V<l>p(f) ■ V<l>^_p(rO (4.12) 
For these equations to be vahd, one requires 

|V$„+i(r-0| « |V$„(r-)| (4.13) 
We retain only the first two terms, and express the lowest-order term in the form 

Uo{r) = exp[$o(r)] (4.14) 
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then equation ()4.9|) is equivalent to the free-space wave equation 

Auo{f) + kluoir) = (4.15) 

Now we set 

W^(r] 

$i(r-) = ^ (4.16) 

By inserting this into equation ()4.10|) and by applying ()4.15|) . we need to obtain 
the solution of 

AWi{f) + klWi{f) = -2klni{f)uo{f) (4.17) 

The solution to this linear differential equation is then obtained with the Green's 
function for free-space propagation as 

$i(f) = / dV'uo(f')nAf')G(\f- ^1) (4.18) 

Uo{r) J 

where the integration is over the source volume denoted by primed coordinates. 
We have obtained the Green's function in Chapter 1, which is 



= (4.19) 



Thus the perturbed field is 



= / dV'u,{^)n,{f^ f''^f_'^\. ""^'^ (4.20) 
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By ()4.8p . the solution of ()4.4|) to the 2nd order is then 

u{r) = Uo{r)exp[5^i{r)] (4.21) 

where $i(r) is given by ()4.20|) . 
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Chapter 5 

Rytov Approximation II - Green's 
Function of Helmholtz Equation 

5.1 Green's Function of Helmholtz Equation by Rytov Ap- 
proximation 

Following the previous chapter, the Helmholtz equation in turbulence for one com- 
ponent is 

Au{f) + kln^{r)u{r) = (5.1) 

We also assume that the magnitude of the inhomogeneity is small, so that the 
refractive index is 

n{f) = I + 6ni{f) (5.2) 
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where Sni{r) << 1, and 6 was inserted into this expression to show smallness. 

Then the Green's function G{f,^) for ()5.1|) satisfies the following equation 

AG(f,0 + fcoV(r)G(f,e) = -5(|f - |j) (5.3) 
In the Rytov method, the solution is expressed as 

G(f,0 = exp($(r,0) (5.4) 
which leads to the following equation 

e*('''«)[A$(f, e) + V$(f, ■ V$(f, 0] + A;2n2(r)e*("^'«) = -6i\f- ej) (5.5) 
or equivalently 

A<l>(f, + V$(f, ■ V$(f, = -fco^^(rO - e^*("'«^5(|f - ^1) (5.6) 

As in chapter 2, we assume that the solution of ()5.(jj) can be written as a power 
series in 6, 

$(f, e, (5) = $o(r, 6 + 5$i(f, + 5'$2(r, + • • • (5.7) 
Notice that the refractive index is given by ()5.2|1 . then the right hand side of 
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equation ()5.6|) could be written as 

i?i/5 = -A;^(l + 25ni(r0+5V(^)-e"*°^'''^'^(l-5$i(r,O + O(5'))5(|r-el) (5.8) 

where we have approximated $(r, ^) by the first two terms in ()5.7|) . and apphed 
Tayler's expansion on e"^*^^'"'^. 

By inserting ()5.7|1 into ()5.6|) and by separating the equations based on the power 
of 6, we obtained the equation for 6^ as 

A$o(r,0 + V$o(r,3 ■ V<l>o(f,e) = -kl - e-*°(^^'«^5(|f - C|) (5.9) 
Now we set 

Go(r-,0 = exp($o(r,6) (5.10) 
then equation ()5.9|) is equivalent to 

AGo(f,0 + ^oGo(r,0 = -Si\r- ?!) (5.11) 

which is the equation to solve for the Green's function for the free-space Helmholtz 
equation. 
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It has been showed in chapterl that 

G'o(r,0 = -— -1 (5.12) 

47r|r - ,^1 

The equation for 6^ is 

A$i(f,0 + 2V$o(r,0-V$i(f,0 = -2A;gni(f) + e-*«(^^'«^)<l'i(f,0<^(|r-e1) (5.13) 
If we set 

W,{f, = $i(r-', OGoir, (5.14) 

then 

AW,{f, = Go{f, A^-ilf, + 2Go{r, OV$o(r, " V$i(r-', + $i(r, O^Go{f, 

(5.15) 

Now by combining (j5.11|) . (j5.13p and (j5.15|) . we obtain 

AW,{f,0 + klW^{f,0 = Go(f,0(A<l'i(f,0 + 2V<l'o(r-',0- V$i(f,0) 

+<l>i(f, Oi^Goir, + k'oG,{f, 0) 

= Go{f,Oi-^klm{r) + e-*»(^^'«)$i(f,e>(|r-'- ^1)) 
-$i(r,e>(|r-el) 

= -2Go(f,e>oV(r1 (5.16) 
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The solution to this hnear differential equation with constant coefficients is ob- 
tained with the Green's function for free space propagation as 



1 



2Go{(:,OKn,{OGo{fX)dC (5.17) 



where the integration is over the source volume denoted by primed coordinates, 
and Go is given by ()5.12p . 
For simplification, we denote 



then the Green's function for (j5.ip to the second order is 



= Go(f,Oe'-^''o"^(^^)^(^^'«^'^^^^^ (5.19) 



5.2 Reciprocity of the Green's Function 

Theorem The Green's function for the 3-D Helmholtz equation is symmetric. 
In other words, if the function G Satisfies the Helmholtz equation 



AG{f, fo) + kln^{r)G{f, fo) = -5{f - fo) 



(5.20) 
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where /cq is the wavenumber, n(r) is the refractive index, and 5{-) is the volume 
impulse function, then 

G{n,fo) = G{fo,ri) (5.21) 

Proof: We shall assume that G satisfies a homogeneous boundary condition at 
infinity, i.e. 

\G{f,fo)\^0 as |f-fo|^oo (5.22) 
From equation ()5.20|) we have that 

G{f,fi)AG{f,fo) — G{f,fo)AG{f,fi) = —[6{f—fo)G{f,fi) — 6{f — fi)G{f,fo)] 

(5.23) 

Let S denote the sphere of infinite radius centered at tq, and V denote the volume 
enclosed by S, then we may integrate ()5.23p over V and obtain 

J^[G{f, n)AG{f, fi) - G(f, fo)AG(f, n)]dV = G(fi, fo) - G(fo, n) (5.24) 
Applying Green's theorem on the LHS of (j5.24j) . we obtain 

J^[G{f, n)VG{f, fo) - G{f, fo)VG{f, n)] ■ dA = G{n,ro) - G{fo, n) (5.25) 
and from the homogeneous boundary condition (j5.22p . we obtain 



G{fi,fo) = G{fo,fi) 



(5.26) 
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End of proof. 
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It's easy to check that the Green's function ()5.19|) we have obtained by Rytov 
method satisfies symmetry property, since 



gifco|r-C| 



G'o(r,0 = -— (5.27) 



is a symmetric function, thus K{f,^,Q = K{^,f,(), where K is defined in (I5.18|) 
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Chapter 6 

Parabolic Approximation 

6.1 Solution of Helmholtz Equation by Parabolic Approx- 
imation 

As in the previous two chapters, the starting point of all propagation theory is the 
Helmholtz equation 



where kg is the wave number, and n(-) is the refractive index, for which we can 
assume the form 



Au{r) + kgn^ {r)u{r) = 



(6.1) 



n(r) — 1 + Sni{r) 



(6.2) 



Next we choose a propagation direction: say, X-axis, and write 



f — ix + jy + kz 



(6.3) 
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p = jy + kz (6.4) 

where i, j, k are orthonormal unit vectors in 3-D space, and x, y, z are scalars. 
Then we let 

u{x,p) = v{x,p)e'-''^ (6.5) 
and substituting in the Helmholtz equation 1)6.11) we obtain 

Fp'u (~)v 

— + 2ik^— + kl{n^ - l)v + ApV = Q (6.6) 



where A^^; = {9^ + q^)v- 



Equation ()6.6|) is equivalent to 



d'^v , dv d ,dv , , , 
— + 2ikoj- = j-{j- + 2ikov) 6.7 
ax^ ox ox ox 



We assume that v is slow-varying in x, so that 



dv 

— + 2ikov ^ 2ikov (6.8) 
ox 



or that the |^ term can be neglected. Then we have the parabolic approximation 



2ikQ— + kl{n^ - l)v + Ap^; = (6.9) 
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In Rytov's method the solution is expressed as 

v{x,p) = exp[^{x, p)] (6.10) 

This leads to the following equation 

^iko^ + Ap$ + Vp$ ■ Vp$ = -klin^ - 1) (6.11) 

where V,$ = (f,f). 

In the perturbation method, we assume that the solution can be written as a power 
series in 6, i.e. 



$(x, p, 6) = $o(x, p) + 5<l>i{x, p) + 52$2(x, p) + (6.12) 

By inserting (j6.2p and (j6.12p into (j6.11|) . and by separating the equations based 
on the power of 5, we obtain the equation for 5^ as 

2iA;o^ + A,$o + Vp$o ■ Vp$o = (6.13) 

Now we set 

^0(2:, p) = e^^^"'") (6.14) 



then equation ()6.13p is equivalent to 
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22A;o^ + Apvo = (6.15) 



The equation for is 



2^A:o^ + A^<l>i + 2V,<l>i • = -2n,kl (6.16) 



If I set 

Wiix, p) = $i(x, p)vo{x, p) (6.17) 
Now by combining ()6.15j) . ()6.16|) and (j6.17j) . we have 



^W^ 
d 

2iko^vo + 22A;o$i^ + voAp^i + 2Vp$i • VpVo + ^iA^dq 
-2klniVo (6.18) 



The solution to this equation is obtained with the Green's function for free space 
propagation as 

vo{x, p) 

° J vo{a,e)ni{a,e)Gp{x,p;a,e)dad9 (6.19) 



2/un 



Voix,p) 
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where the integration is over the source volume, and Gp is the Green's function for 

dnmi). 

The Green's function has been obtained by Ishimaru [21] as 



Gpix,p;a,e) = (6.20) 



Thus the perturbed field is 



<l>i(x,p) = — ^ / v^{a,e)nr{a,e)— -dadO (6.21) 

fo(x, p) J Ati[x — a) 



The solution to the 2nd order is then 

v{x,p) =VQ{x,p)e^'^'^''^P^ (6.22) 

and 

u{x, p) = vix, p)e'^" = Vo{x, p)e'5*i(-''')+^'=^ (6.23) 

We specialize next to two cases. 
Case 1. Plane Wave 
Here we take 

uoix,p) = e''^ (6.24) 

and hence 

vo{x,p) = l (6.25) 
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since Uq = VqC^^^ . Then equation ()6.2H) yields 



/g 2(x-a) 
ni{a, 0)-— -dadO 
4:7r[x — a) 

kl fni{a,e) jh^ 



2n 



I'h^e'^dade (6.26) 
J X — a 



Case 2. Beam Wave 
For the beam wave case 



«o(0,p) = exp(-^) (6.27) 



and hence by Ishimaru [21], we have 



Uo{x.p) = j^. — exp{-—^^. — -) (6.28) 
1 + ixa 2(1 + ixa) 



which is equivalent to 



, , 1 , akp^ , ,„ 

Vo[x,p) = — — — exp{- . — -) 6.29 
1 + txa 2(1 + txa) 



Hence 



ifc|p-9r 

(Pi[x,p) = — ; — / expi z -)ni[a,0) — ; -dadU 



k?. . , , , akp'^ 
— (1 + ixa)exp[- 



27r' ' '2(1 + ixa)' 

J {l + iaa){x — a) 2(1 + zaa) 2{x — a) 

(6.30) 
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6.2 Green's Function by Parabolic Approximation 

The Green's function p; a, 6) for ()6.9|1 satisfies the following equation 
d 

2iko—G'{x, p; a, 9) + kl{n^ - l)G\x, p; a, 9) + ApG'{x, p; a, 9) = -S{x, p; a, 9) 

(6.31) 

where 5{-] ■) is the Dirac delta function. 

In Rytov's method, the solution is expressed as 

G"(x, p; a, 9) = exp{<^{x, p; a, 9)) (6.32) 

which leads to the following equation 

'2iko^ + Ap$ + Vp$ • Vp$ = -kl{n^ - 1) - e-^6{x, p; a, 9) (6.33) 

We assume that the solution to ()6.33p can be written as a power series in 5, 

<l>(x, p- cT, 9- 5) = %{x, p; a, 9) + 5$i(x, p; a, 9) + 5'^^2{.x, p;(r,9) + ... (6.34) 

Notice that the refractive index assumes the form n{-) = 1 + 6ni{-), the right hand 
side of equation ()6.33p could be written as 

RHS = -kl{25ni + d^nj) - e-*°(l - + 0{5^))5{x, p; a, 9) (6.35) 
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where we have approximated $ by the first two terms in ()6.34|) . and apphed Tayler's 
approximation on exp{—S^i). By inserting ()6.34|) into ()6.33p and by separating 
the equations based on the power of 6, we obtained the equation for as 



+ Ap*o + Vp$o ■ Vp$o = -e-*°(^(x, p; a, 9) (6.36) 



Now we set 

Go(x, p; a, 6) = exp{^oix, p; a, 6)) (6.37) 



Then equation ()6.36p is equivalent to 



d 

2iko—Go{x, p; a, 9) + ApGo{x, p; a, 9) = -6{x, p; a, 9) (6.38) 



As we know, the Green's function for free-space propagation has been obtained as 



Go{x,p-a,e) = / (6.39) 
4:T{[x — a) 



The equation for 5^ is 



2«A;o^ + A^$i + 2Vp$i ■ Vp$o = -2A;>i + $ie-*"5(x, p; a, 9) (6.40) 



If I set 

Wi{x,p]a,9) = (^i{x,p;a,9)Go{x,p;a,9) (6.41) 
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dx 



2iko^ + ApWi 



-2klniGo (6.42) 



where we have apphed ()6.38p and ()6.40p . The solution to this hnear differential 
equation is obtained with the Green's function for free space propagation as 



Go{x,p; a, 6) 
Go{x,p; a, 0) 



Go{t, 7; a, 6')ni(r, -f)Go{x, p; r, -f)dTd-f 

(6.43) 



where the integration is over the source volume, and Go is given by ()6.39j) . 
For simplification, we denote 



7^// n ^ GoiT,-f;a,0)Go{x,p;T,j) 

K {x,p;(T,e;T,-f) = — (6.44) 

Go{x, p; a, 9) 



then the Green's function for (j6.31|) is 



G'{x,p;a,9) = exp{^{x, p;a,6)) 
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= Go{x, p;a,9)exp{S J 2klni{T,'y)K'{x, p;a,6;T,^)drd'y) 

(6.45) 

It's easy to see that the Green's function ()6.45|) we have obtained by using parabohc 
approximation satisfies the symmetry property. 
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Chapter 7 

Time Reversal Application 

An important phenomenon in wave propagation in non-homogeneous media is re- 
lated to time reversal of the wave field. In the time reversal procedure, the wave 
received by an antenna (receiver-emitter) array is recorded and then re-emitted 
into the medium time reversed, that is , the tails of the recorded signals are sent 
first. The time-reversal procedure is equivalent to phase conjugation on the spatial 
component of the time-harmonic wave field. 

A cartoon of a time reversal experiment is given in Figure ()7.1|) 
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Minor 




|-« z >| 

Figure 7.1: The time reversal procedure. A pulse is emitted from a source with central 
wave length Aq. The transmitted field is recorded, stored and time reversed at the mirror 
of size a at z distance away, and then sent back toward the source. It refocuses on the 
spot of size r*. 

To simplify the problem, let the phase-conjugated mirror be located at z = 0, 
and the source at the parallel plane 2;-distance away. The aperture function of the 
mirror is assumed to be the indicator function xa, where the set A represents the 
physical boundary of the mirror. 

Let Gniro, ri) be the Green's function for the 3-D Helmholtz equation in random 
media, then G{fo,fi) defined in (I5.19|) is an approximation. 
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In section 3.2, we stated and proved that G{fo, fi) satisfies the symmetry property 

Gifo,n) = G{n,fo) (7.1) 
The wave field received at the mirror is given by the super-position integral 

where the integration is over the source volume denoted by primed coordinates, 
and \l/o(^s) is the input field. 



After phase conjugation and back-propagation we have at the source plane the 
wave field. 

^^(r) = J G{f, fm)G{f^, rs)xA{rm)M^df^dfs (7.3) 
By inserting ()5.19|) into ()7.3|) . we obtain 

■ X Air m)'^ air's) dfmdfs (7.4) 
where Goif,0 = and miX) = ^^^flf^- 
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Let E{-) denote the ensemble average, we can evaluate the field received at the 
source plane by 



£;(^^(f, ni)) = E{ / Go{f, f,n)Go{f^, ^^)e5/ 2/c2„i(a(i^(r,r„,C)+X(r^,r- ,a)dC 



■XA{rm)'^o{rs)df„^dfs) (7.5) 



Go{f, fm)Go{fm, f,)E(e^/ 2^oV{a{i^(r>--™,C)+/<(r™,r.,C))<iCj 



■XA{rm)'^o{rs)dfmdfs (7.6) 



Go{r,f^)Go{fm,rs 



■e 



452fc4 J E(ni((i)ni{(2))iK{r,r,n,Ci)+K{rrn,rs,Ci))iK{r,r^ 



■XA{rm)'^o{rs)df„^dfs (7.7) 
where E{ni{(i)ni{(2)) is measured by experiments. 

Experimentally, it is found that when the wave amplitude \Ef^ is recorded, either 
photographically or by striking an observation screen, we record only the intensity 
/ = Hence, from ()7.4|) we have 

/(^ = |$^(f)|' 



drmidfsi / Goir,f^2)Goir^2,rs2)xAirm2)'^oirs2) 
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= j Go(r, rmi)Go{fmi,fsi)Go{f, rv„2)Go(^m2, r;2)xA(rv„i)XA(rVn2)*o(^si)*o(^s2) 

J 2k^ni{C)[K{r,rrnl,0+K(frnl,rsl,0+K{r,rrn2,0+K{rrn2 

(7.8) 

where Go(f,a = and K(f,e,C) = ^"^ggf^ - 

By evaluating the ensemble average of the intensity, we obtain 

E{I{r,ni)) ^ E{j Go(^ ^mi)Go(f^i, r^^Gol^ ^m2)Go(r^2, 

+K{f, rrn2, + K{fm2, ^2, 0]dQdrm\drsidfm2dfs2) 

= j Go{f,frnl)Go{fral,fsl)Goir,frn2)Go{frn2,rs2)XAir^^^^^^ 

■^^{¥^Mrs2)E{exp{5 j 2/^oV(C)[i^(r, rVni, C) + ^(^ni, r,i, C) 

+K{f, frn2, + K{fm2, ^2, 0]dC}) df mid fsld'rm2dfs2 
= j Go{f,fmi)Go{f„ii,fsl)Go{f,frn2)Go{frn2,rs2)XA{rml)XA{rm^ 

■^^Mrs2)exp{46^k', J E{n,{Ci)ni{C2))[K{f, fmi, Ci) 

+K{fml, fsi, Cl) + K{f, fm2, Cl) + K{fm2, rs2, Cl)][K{f, C2) 
+K{fml, fsi, C2) + K{f, fm2, C2) + K{fm2, rs2, C2)]c^ClC?C2} 
dfmldfsldfm2dfs2 (7.9) 
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where the data we need to evaluate E{ni{(i)ni{(2)) is measured by experiments. 



To illuminate this method, let's consider one of the simplest cases. Assume that 
"^ol') and Xa{') are Dirac delta functions, which means that both the source and 
the phase-conjugating mirror is a single concentration point. Thus we have 

^m(rm) = S{fm)G{f,n, r,) (7.10) 



and 



*^(f,) = Gifs,rm)Gir^,fs) (7.11) 

where fg and fm are position vectors for the source and the mirror respectively. 
By inserting ()5.19j) into ()7.1H) . we obtain 



= |G'o(r,,f^)|2e'^/4'=oV(C>e(i^(r-,r„,C))dC (7.I2) 



where Goir.f) = °J , and K(f,FX) = ^ ■, and we have applied the 

symmetry property of Go- 

We can also obtain the intensity as 



/(f,) = |vi/^(r-;)p 

= vl/^(r-;)^^ 
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= \Go{f,,f„,)\''e^I^''>'^^'^^'^^^^^^^"^'^^^^ (7.13) 
By evaluating the ensemble average, we obtain 
E{I{fs,ni)) = £;(|G'o(f.,f^)|V/«^om(ORe(i^(r./^,C))<iC) 

= |G'o(f.,nn)|'e6^^'-'*^o/^<^(------.0)^dC (7.14) 

where we have taken the Gaussian expectation and assumed that ni(-) is white 
noise, so that £;(ni(Ci)ni(C2)) = (^'^^{Ci - (2)- 



CHAPTER 8. 



3-D Apodization Problem In Turbulence 



50 



Chapter 8 

3-D Apodization Problem In 
Turbulence 

In this chapter, we are concerned with the determination of the distribution of 
light over the exit pupil of an optical system required in order to achieve a desired 
distribution of illuminance over a given plane in the image field. This problem is 
known as the Apodization Problem, and we are particular interested in the deter- 
mination of that amplitude distribution over a circular pupil which maximizes the 
fraction of the total energy that lies in a prescribed circle in the image plane. 

Let X = (x^, x^) be the radius vector in the plane of the exit pupil from the optical 
axis to an arbitrary point in that plane; let ^ = (^^,^^) be the radius vector in 
the image plane from the optical axis to a point in the image plane(Figure (|8.1|) ). 
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I z 

Figure 8.1: A configuration of the Apodization Problem. 

Then by the extended Huygens-Fresnal principle, the light amplitude A{£) in the 
image plane is given by the superposition integral 

m = / Gh(0, X- z, i)T{x)dx (8.1) 

J\x\<a 

where Gh{0,x; z,^) is the atmospheric impulse response (Green's function), T{x) 
is the light amplitude in the circular exit pupil of radius a, and z is the distance 
from the pupil to the image plane. 

Our apodization problem then requires finding the function T{x) for which the 



CHAPTER 8. 3-D Apodization Problem In Turbulence 52 
ratio 

VI-4K)P<i« 

is a maximum. 

Here b is the radius of the circle in which the illuminance is to be maximally con- 
centrated, and the circle is centered on the optical axis. 



We assume the optical wave propagates through a thick slab of turbulence air, 
where Rytov's method could be applied to approximate the Helmholtz equation. 
Prom chapter 4, we have derived an explicit formula for Gh by Rytov's Approxi- 
mation. 



where 



i^(nc:o^ ^°'^ff^' (8.4) 



and 



gifco|r-g| 



^0(^0 = -—^ (8.5) 
47r|r - ,^1 

Denote the top of ()8.2|) by M, then M is a function dependent of random variable 
ni(C), i.e. 

M(ni)= / (8.6) 
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= I GniOjXi] z,£)T(xi)GHiO,X2; z,^)T(x2)dxidx2d^ 

J \£\<b.\x9\<a,\xi \<a 



T{xi)T{x2)Go{0,xi; z,^)Go{0,X2; z,^) 

\^\<b,\x2\<a,\xi\<a 



■e 



5 J 2fc2ni(C)/r[{0,xi),(^,a,C]rfCg5 / 2fc2„i(c)x[(0,a-2),(^,0,C]dC^-^^-^^^ 



T(fi)r(x2)G'o(0, xi, z, 0^0(0, X2; Z, ^) 

\^\<:b,\x2\<a,\xi\<a 



.g5 / 2fc2„i(^--){i^[(o,xi),(^,0,cl+if[(0,X2),(2/),Cl}'iC^f^(^f2rf^ (8.7) 

where we have assumed «.i(C) to be real. 

Let E{-) denote the ensemble average (expectation), we need to evaluate the statis- 
tics of the energy received at the source plane. 

E{M{m)) = E{[ ^ r(f 1)7X^^0(0, xi; z, 0^0(0, X2; z, 



■e 



5 J 2fc2ni(0{X[(0,xi),(2,0,Cl+X[(0,X2),(^/),Cn<iQ^^t/f2fi6 



T{xi)T{x2)Go{0, Xi, z, 0^0(0, X2; z, ^) 

5|<b, |2?2|<a, |a?i|<a 



_^^^5 J 2fc2ni(a{E-[(0,xi),(2,a,Cl+^[(0,a;2),(;z,0,f]}rfC)^£^tiX2(i^ 

r T{x^)TXS2)GoiO,xi,z,C)GoiO,X2;z,Oexp{Ad'k^o 

^ IfKfi.la;'! I<a,la;i \<a 



E{n,{Ci)nim{Km fi), {z, OXi] + K[{0, X2), {z, 0, Ci]} 



■{K[{0, fi), {z, 0, C2] + K[{0, X2), {z, 0-X2]}dCidC2)dx,dx2d^ (8.8) 



CHAPTER 8. 3-D Apodization Problem In Turbulence 54 

Assuming ni(-) is white noise, i.e. i?(nl(Cl)?^l(C2)) = o^^^iCi — we can simplify 
dHSl) as 

E{M{m)) = I ^ r(fi)T(^G'o(0,fi;;2, 0^0(0, fa; ^,0 

l?l:^^il^2 |<a, 1^1 1 <a 



r(f i)r(f2)G'o(0, fi; 2;, ^)G'o(0, ^2; 2, 

g|<6, |x2 |<<i, |xi |<a 

g452<T2fc4 J (i^[(0,a;i),(^/),Cl2+2i^[(0,a;i),(;^/)/]X[(0,X2),(^/),Cl+i^[{0,X2),(z/^ 

T(f 1)^0(0, fi; f)e^^''^''=oV^[(o.-i).(--a.cT«(C 

|4]<ft,|x2|<a,|5i|<a 



.g85V2fc„Vi^[(oA),(^,a,ax[(o,x2),(.,a,c]<ic^^^^^2^l' (s.g) 

Since << 1, we can rewrite the last term inside the integral of ()8.9|) by Tayler's 
expansion as 



= l + S5'aX j K[{<d,x,),{z,0X]K[{Q.X2).{z,0X]dC + O{6') (8.10) 

Denote 
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= ( \ ( D{x,Odx\'^d^ 

J\i\<b J\x\<a 

+ Sd'aXf. \f D{x,OK[iO,x),iz,0,C\dx\-'dCdC (8.12) 

J\^\<b J\x\<a 

If the air-density inhomogeneity is extremely small, which means S'^k^ << 1, then 
the second term is negligible compared to the first one. Consequently we have 

E{M{ni))= [ \[ D{x,i)dx\'^di (8.13) 

J\i\<b J\x\<a 

This result can be obtained equivalently by taking the expectation of the light 
amplitude A{S,) first, and then by evaluating the energy concentrated in the image 
circle from the expected light amplitude E{A{S,)), since 

= E{[ GH{0,x;z,i)T{x)dx) 

J \x\<a 

= E{f GoiO, x; z, Oe'^/2'=o"^«^)^[(°'"^)'(^'«^'^^'^^T(f)(if) 
= / r(x)G'o(0, x; z, f)E(e^/2A^g"i(C)^[(o.-).(-'a,C]<iC)rf^ 

J |x|<a 

= / r(x)G'o(0, X; Z, ^)e^^^''o I E{ni(a)ni{C2))Kl{0,x)Xz£),Ci]Kl(0,x)Xz,(),C2]dC^^^^ 
J\x\<a ' ' ' 
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x\<a 

D{x,^)dx (8.14) 
where 

D{x, = r(x)Go(0, f; z, ^)e^s'-'''o I Klio,m-,iU?d( (g 
which is the same as (jS.llll . 

Therefore the energy received in the image circle is 



M' = ( \E{A{0)\Hi 

J\i\<b 

= I \ I D{x,^)dx\^di (8.16) 

J\£\<b J\x\<a 



which is identical to ()8.13|) . 

Furthermore we can assume that the energy is conserved through wave propaga- 
tion. Therefore by similar calculation, we obtain the following equation for the the 
bottom of (HIS) 

/ \A{i)\'di= I I / D{x,i)dx\'di (8.17) 

JR? JR? J\x\<a 

Thus our original apodization problem is reduced to the one that requires finding 
the function T{x) for which the ratio 



I !\x\<aD{x,i)dx\'^d^ 
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is the maximum, where D{x,0 := T{x)Go{0,x; z,i)e'^^'''''''^ I and 

If A is of total energy E, since the energy is concerved, we have 

[ I / D(x,e>dx\^di^ I \Tix)\'^dx (8.19) 

Jb? J\x\<a J\x\<a 

whereas the energy in the image circle is 



D{x,i)dx\^di 

|^|<6 ■ J\x\<a 



D{x,^D{y,Odxdydi 

|^|<6 ^|a;|<a J\y\<:a 



/ / K,{x,y)T{x)T{y)dxdy (8.20) 

J\x\<a J\y\<a 



where 



Ks{x, y)= [ Go{0, X- z, OGo(0, y; z, f)e4^^-^'=oV(^[(o--l-(--0-CT+^[(o,5)-{^>0>ClVQ^-' 

(8.21) 

Our apodization problem therefore requires finding the function T(x) for which 
the ratio 



_ I\x\<al\y\<a^s{x, 

is the maximum. This maximum is Aq, the largest eigenvalue of the integral equa- 
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X'^{x)= f K,{x,y)^{y)dy, \x\ < a (8.23) 

J\y\<a 

with kernel 

K,{x, y)= I Go(0, X- z, 0^0(0, ^, Oe''''^''» /(K[(o,^),(.,a,cT+i^[(o,j;),(.,a,ClVcrf^-' 
J\i\<b 

(8.24) 

Somewhat simpler than ()8.23p is the integral equation 

a^{x)= [ Go{0,x;z,Tf}e^^''''''olKKo,xUz,rf),^^dc^^^-^^^ (g25) 

J\rf\<a 

Whe shall show below that the solution of this equation is completely equivalent 
to the solution of (jHI2Sl)- 

From the symmetry of the input domain, it follows that if \l/(x) is a solution 
of (IHI2S1), so also is ^(-x), so that both ^'e(^) = ^(x) + ^(-f) and ^o(^) = 
\I'(a;) — \E'(— x) are solutions as well. The eigenfunctions of (j8.25|) can be chosen to 
be either even or odd functions of x. 
The complex conjugate of ()8.25p is 

J \f[\<a 

Multiply (jg^ by \l/(f) and integrate over |f | < a. Multiply (|g^ by \l/(f) and 
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integrate over \x\ < a. Combining these equations, we find on using the symmetry 
of the integral domain that 

(a±a) / \I/(f)\l>(f)(if 

•'\x\ <a 

(8.27) 

If then \l/ is even, by choosing the negative sign in this equation, one obtains 
q; — a = 0, whereas if ^ is odd, by choosing the plus sign, one finds a + oi = 0. The 
eigenvalues of ()8.25p associated with even eigenfunctions are real; the eigenvalues 
of ()8.25|) associated with odd eigenfunctions are pure imaginary. If follows then 
that (j8.25|) is equivalent to the pair of equations 

J\x\<a 

(3o^,{x)= I /m(G'o(0,x;z,r/)e^'^''^''=o/^[(o.-^-(-.'?)-CTrfC)^^(^)rf^ (8.29) 

J\x\ <a 

in which f3e and /So are real. These equation have real symmetric kernels and we 
can fall back on the extensive theory in the literature treating such equations. 
It follows then from [22] that the eigenfunctions of (j8.25|) can be chosen real, 
orthogonal and complete in the class of functions square-integrable in the input 
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J\fj\ <a 



\v\<a 



y\<a 



/ Kr.{x,y)<il{y)dy (8.30) 

J\y\<a 

A'=|aP (8.31) 



where 



J\fJ\<a 

(8.32) 

which is ()8.24|) in shghtly ahered notation when a = b, since Go{0,fi;z,y) = 
Go{0, y; z, ff). When a b, it's a scaled version of (j8.24p by paraxial approximation. 
Since the solution of ()8.25|) are complete, it follows that they are also a complete 
set of solutions of ()8.30|) . As we asserted, to solve ()8.3Up . it suffices to solve ()8.25|) . 
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Chapter 9 

Conclusion and Discussion 

In the first part of this thesis, we have been deahng with the 3-D Helmholtz equa- 
tion which is fundamental for all propagation theory. The Green's function for 
free-space propagation could be easily solved with specified boundary condition 
(Sommerfeld radiation condition). 

In a turbulent medium where the magnitude of the air-density inhomogeneity is 
small, Rytov's method could be applied to approximate the Helmholtz equation 
for the frozen atmosphere. Within a limited but useful region of validity, both the 
solution and the Green's function of the Helmholtz equation are obtained explicitly. 

Furthermore, by using the parabolic approximation, a simpler form could be ob- 
tained for both the perturbed wave field and the Green's function. In two spe- 
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cialized cases for the initial condition which are the plane wave case and beam 
wave case, the solution could be expressed in a form that numerical computation 
is easier to conducted. 

In the second part of this thesis, through direct application of the extended Huygens- 
Fresnel principle, a general expression has been derived for the wave field received 
in the object plane for two optical problems - time reversal and apodization prob- 
lem. In order to evaluate the ensemble average of the field and the intensity, we 
assume the refractive index disturbance to be white noise and take the Gaussian 
expectation on the inhomogeneous term. By doing so we end up with superposi- 
tion integrals which can be computed by using numerical integration techniques. 

At the very end of Chapter 8, we obtain an energy ratio (|8.22|) which is very similar 
to the one for the free-space apodization problem except for an exponent term in- 
volved with the air-density inhomogeneity. The apodization problem for free-space 
propagation, which means 6 = in our case, has been well studied t^^ll^^H^^l, and 
the free-space input and output eigenfunctions are prolate spheroidal wavefunc- 
tions. Therefore an immediate research problem based on this thesis is to build a 
connection between our results through Rytov's method and previous analysis of 
the analytic solution of the apodization problems. 
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To obtain the numerical results for these optical applications, we finally need to 
do integrations in 2-D or 3-D domains. In doing so, the parabolic approximation 
is helpful since it can separate the first dimension with the other two dimensions 
and thus reduce the complexity of the integrals. We can also specify the initial 
conditions in particular cases such as the plane wave and the beam wave. 
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